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ABSTRACT 

An efficient approach for simultaneous aerodynamic analysis and design optimization is 
presented. This approach does not require the performance of many flow analyses at each 
design optimization step, which can be an expensive procedure. Thus, this approach brings 
us one step closer to meeting the challenge of incorporating computational fluid dynamic codes 
into gradient-based optimization techniques for aerodynamic design. An adjoint- variable method 
is introduced to nullify the effect of the increased number of design variables in the problem 
formulation. The method has been successfully tested on one-dimensional nozzle flow problems, 
including a sample problem with a normal shock. Implementations of the above algorithm are 
also presented that incorporate Newton iterations to secure a high-quality flow solution at the 
end of the design process. Implementations with iterative flow solvers are possible and will be 
required for large, multidimensional flow problems. 
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Adjoint variables 
Local cross-sectional area 
Design variables 
Specific internal energy 
Specific total energy 

Quasi-one-dimensional Ruler equation flux vector 

Number of constraints 

Pressure 

State variables (aerodynamic field variables) 

Nonlinear state equations or residual of flow equations 
Penalty function 

Quasi-one-dimensional Euler equation source vector 
time 

Flow speed 

Nozzle axial coordinate 

Step size (during one-dimensional search) 

Ratio of specific heals 
Change in quantity 
Merit function 
Density 

Side-constraint functions 
Objective function 


Subscripts 

i ith side constraint 

j x-coordinate discretization index 

o Objective function 

Superscripts 


n, m 


Subiteration index 
Iteration indices 



T 


Matrix transpose 
Target values 
Approximate operator 
New (or updated) quantity 

Abbreviations 

ADS Automated Design Synthesis 

CFD Computational fluid dynamics 

CPU Central processing unit 

FD Finite difference 

NR Newton-Raphson iterations 

QA Quasi-analytical 

SAADO Simultaneous aerodynamic analysis and design optimization 

2.0 INTRODUCTION 

A typical aerodynamic problem can be highly nonlinear, and its solution can consume a great 
deal of computer resources in terms of both memory and central processing unit (CPU) time. 
This demand poses a challenge for incorporating advanced computational fluid dynamics (CFD) 
software into gradient-based optimization techniques for aerodynamic design. In this study, 
a strategy is proposed whereby the objective function and constraint violations are reduced; 
simultaneously, the aerodynamic field equations are solved (i.e., the final steady-state numerical 
solution of these equations is achieved). With this approach, the aerodynamic field variables 
are treated as additional design variables, and the aerodynamic residual equations are considered 
to be additional equality constraints. As a consequence, the steady-state aerodynamic equations 
are satisfied only at the final optimum solution. With this new approach, overall computational 
expense is expected to be reduced in comparison with the conventional optimization method (for 
which many repeated steady-stale aerodynamic analyses are performed during the optimization 
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process). Two quasi-one-dimensional nozzle-flow problems are considered to demonstrate the 
initial numerical testing of the scheme. In addition, the feasibility of extending the methodology 
to multidimensional aerodynamic problems is demonstrated. 

Other recent studies involve the development of strategies for simultaneous aerodynamic 
analysis and design optimization. For example, the algorithm developed by Rizk [1, 2] first 
updated the design variables at each design optimization iteration with a gradient-projection 
method and then updated the field (state) variables with several flow-solver iterations. The 
gradients of the objective and side-constraint functions were calculated for the current values 
of the flow variables by the method of finite differences. The inverse design scheme of Drela 
[3] employed a streamline formulation of the Euler equations; in addition, the locations of the 
streamlines that define the shape of the airfoil were taken as the design variables. The optimum 
design and the corresponding flow solution were simultaneously obtained by solving a global 
system of equations that included (at each optimization iteration) the linearized flow equations 
and the necessary conditions of the minimization problem. 

Design optimization techniques are generally iterative in nature; the fact that simultaneous 
analysis and design is efficient only when the method of analysis is also iterative is noted by 
Haflka [4, 5]. Recent works by Orozco and Ghattas [6-8] and Ta’asan et al. [9, 10] investigated 
techniques to improve the computational efficiency and storage requirements for simultaneous 
aerodynamic analysis and design optimization. The former method examined the sparsity of the 
Jacobian and Hessian matrices; the latter method took advantage of the multigrid method. 
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Two important concerns exist for the present study: the increased size of the problem 
(because the aerodynamic field variables are treated as design variables) and the quality of 
the aerodynamic solution at the end of the design process (measured by the h norm of the 
aerodynamic residual equations). The first concern is addressed by employing the adjoint- 
variable method; the resulting equations are similar to those derived in reference [6], The 
second concern is addressed by adding flow-solver iterations between the design iterations in 

the optimization procedure. 

3.0 OPTIMIZATION METHODOLOGY 

A typical engineering design optimization problem attempts to find the best design b that 
minimizes the objective function, satisfies the constraints, and satisfies a set of equilibrium 
conditions (i.e., a set of state equations). Mathematically, this optimum design can be stated as 

"j" *o(Q(b),b) (1) 

which is subject to 

#i(Q(b),b)<0 (i = l,2,...,m) (2) 

where the vectors Q and b are the aerodynamic field and design variables, respectively, and 
and are the objective and side-constraint functions, respectively. In an aerodynamic problem, 
the field variables usually represent the density, velocity components, and pressure, for example. 
The field variables describe the physical behavior of the system to be designed and are related 
to the design variables through the nonlinear state equations, which are defined symbolically as 

R(Q(b), b) = 0 (3) 
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As is characteristic of most realistic aerodynamic problems, the dimension of the vector b is 
of the order of tens to several hundreds; in contrast, the dimensions of the vectors Q and R 
are of the order of thousands to even millions. A solution Q for Eq. (3) for a given b can 
be not only difficult to find but computationally intensive as well. However, possibly hundreds 
of such analyses (i.e., repeated solutions of Eq. (3)) are required for the conventional design 
optimization process. The conventional design optimization process is explained in greater 
detail in the following section. 

3.1 Conventional Approach For Aerodynamic Design Optimization 

First, to gain a better understanding of the new algorithm, the general procedure for the 
conventional approach to a typical aerodynamic design optimization problem is outlined in detail. 
From an initial design b, a typical gradient-based design optimization scheme first solves for 
the field variables Q with the state equations (Eq. (3)) and then evaluates the change in design 
Ab that reduces the change in the objective function and corrects the constraint violations. This 
most favorable Ab is obtained by solving a linearized optimization problem 

A b 'MQ< b ).b)+( lnr )Ab (4) 

which is subject to 

#i(Q(b),b) + Ab <0 (i = l,2,...,m) (5) 

where Eqs. (4) and (5) are linearized approximations of Eqs. (1) and (2), respectively. Note 
that for simplicity only the first-order terms appear in Eqs. (4) and (5). The derivatives of the 
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objective function can be evaluated either in terms of the derivatives ^ or in terms of 
the adjoint variables Ao as 

dtf 0 _ 3*0 dQ d* 0 
db dQ db + db 


or 

d*o _ A T dR 

db ~ ° db + db 

where the derivatives ^ are the solution of the linear sensitivity equation 


_ iJQ = ^5; (8) 

db “ rib 

This linear sensitivity equation results from the differentiation of Eq. (3) with respect to b, and 
the adjoint variables A 0 are the solution of the linear equation 



( 9 ) 


The derivatives of the side-constraint functions ^ can be evaluated with expressions that are 
similar to Eqs. (6) through (9). 

After the most favorable direction of change Ab is found by solving Eqs. (4) and (5) with a 
proper mathematical linear programming technique, the next procedure is to determine the step 
size o, which defines the newly updated design b* as 


b* = b + oAb 


( 10 ) 


'pjg step size o serves as a relaxation factor in an iterative scheme for determining b , it is 
determined by solving a nonlinear optimum design problem (similar in form to Eqs. OH 3 )) 
in which o is the only design variable. The procedure for determining o is referred to as the 
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“one-dimensional search” in many design optimization textbooks. The updated values of the 
field variables Q*(b*) for the new design b* in this procedure is obtained by solving the state 
equations (one solution for each trial value of o); that is, 

R(Q*(b*),b*) = 0 (11) 


Finally, the optimum design is achieved if the Kuhn-Tucker conditions are satisfied; if not, the 
entire procedure (Eqs. (4) through (11)) is repeated until convergence is achieved. 

Generally speaking, the conventional design optimization approach of Eqs. (4) through (11) 
requires on the order of possibly hundreds of design iterations to converge to the solution of 
the original design optimization problem (defined previously by Eqs. (1) through (3)). The 
repeated solution (to steady state) of Eq. (1 1) is computationally intensive, particularly for 
three-dimensional nonlinear aerodynamic problems. This computational intensity has motivated 
several studies (including the present one) with the goal of reducing this computational effort 

[i-io]. 

3.2 Proposed Approach For Simultaneous Aerodynamic Analysis 
and Design Optimization (SAADO) 

3.2.1 Derivation of SAADO Equations 


The proposed new approach reformulates the design optimization problem presented in Eqs. 
(1) through (3) as 


which is subject to 


min 

Q,b 


*o(Q,b) 


( 12 ) 


*i(Q,b)<0 (i= l,2,...,m) (13) 
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and the additional equality constraints 


R(Q,b) = 0 


(14) 


This formulation is referred to as simultaneous analysis and design optimization (SAADO), 
which (in contrast to the conventional approach) now treats the state (field) variables as part of 
the set of independent design variables and considers the state equations to be part of the set of 
side constraints. Because satisfaction of the equality constraints of Eq. (14) is required only at 
the final optimum solution, the steady-state aerodynamic field equations are not solved at every 
design optimization iteration. As a result, the excessively large computational burden of the 
conventional approach may be reduced significantly. However, this advantage is likely to be 
offset by the large increase in the number of design variables and equality constraint functions 
unless some remedial procedure is adopted. 

The new method begins with a linearized design optimization problem (similar in principle 
to that of Eqs. (4) and (5)), which is solved for the most favorable change in the design variables 
Ab, as well as in the state (field) variables AQ; that is 


min 

AQ, Ab 


;iil; 

*o(Q, b ) + -^AQ + 


M 0 

db 


Ab 


(15) 


which is subject to 


* i( Q,b) + ^AQ+^Ab<0 


(i = 1,2, ...,m) (16) 


and 


R (Q, b )+^AQ + fAb = 0 


(17) 
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Equations (15), (16), and (17) are linearized approximations of Eqs. (12), (13), and (14), 
respectively. Note in this formulation that R(Q,b), which is the residual of the nonlinear 
aerodynamic field equations, is not required to be zero until the final optimum design is 
achieved. In comparison with Eqs. (4) and (5), the mathematical linear programming problem of 
Eqs. (15)-(17) becomes more computationally difficult to solve because of the dramatic increase 
in the effective number of design variables and equality-constraint equations. This difficulty 
is overcome by the introduction of the adjoint-variable method to remove AQ and Eq. (17) 
altogether from this linear programming problem. 


The derivation of this technique begins with the introduction of an arbitrary vector A 0 to 
augment the linearized objective function of Eq. (15) by using Eq. (17); the resulting enhanced 
linearized objective function is expressed as 


*o(Q.b) + ^?AQ+^Ab + Aj' 


OQ 


db 


R(Q,b) + ^AQ + fA b 


(18) 


If rearranged, then result is 


♦«(Q,b) + 




+ A 1 — 
f A " <)Q 


) 


AQ + A^'R(Q,b) + 


( 


d * 0 

db 


T dR' 


TA^jAb (19) 


Then, the adjoint-variable vector A 0 can be specified such that the coefficient of AQ vanishes, 
which completely eliminates AQ from the enhanced objective function of Eq. (19); this result, 
of course, dramatically reduces the size of the linear programming problem to be solved. Thus, 
the adjoint equation for A 0 is 



( 20 ) 
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and the resulting equation for the enhanced, linearized objective function becomes 


*o(Q,b) + AjR(Q,b)+ (^ + Aj^Ab (21) 

A similar procedure must be applied to reformulate the linearized side-constraint functions of 
Eq. (16). The final result is 

vpi(Q,b) + Ai r R(Q,b)+^^ + A^'^Ab<0 (i = 1,2, m) (22) 


where Aj are the solutions to the adjoint equations 

1 ' 




(i = 1,2, ..., m) 


(23) 


With AQ and Eq. (17) eliminated, the linearized optimum design problem of Eqs. (15)-(17) 


is reformulated as 


min 

Ab 


*MQ» b) + Aj R(Q, b) + 


( + A T — ^ 

V db + ° db) 


Ab 


(24) 


which is subject to 

*i(Q, b) + A', l 'R(Q, b)+ + A ^) Ab < 0 0 = 1, 2, m) (25) 

After Eqs. (24) and (25) have been solved for Ab (with a suitable mathematical programming 
technique), Eq. (17) can be solved to obtain AQ. A one-dimensional search is performed to 
find the updated values of Q* and b*; the search procedure must solve a nonlinear optimization 
problem of the form 

n,in *o(Cr,b*) (26) 
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which is subject to 


*i(Q*,b + ) <0 


(i = 1 , 2 ,..., m) 


(27) 


and 


R(Q*,b*) = 0 


(28) 


where the step size a is the only design variable. Again, note that the equality constraints (Eq. 

(28) ) are not required to be zero until the final optimum design; violations of these equality 
constraints should simply be progressively reduced during the SAADO procedure. Therefore, 
the updated Q* in this study is defined as Q* = Q + AQ*, which satisfies the first-order 
approximation of Eq. (28) as 

R(Q ’ b) + AQ * + k (aAb) = 0 (29) 

The update for b* is given previously in Eq. (10). Note that in the one-dimensional search, Eq. 

(29) need to be solved for each trial value of a; alternatively, a one-time solution (per design 
cycle) of the following two equations is performed: 


<7R A _ 

«q AQi 


JR 

<)b 


Ab 


(30) 


and 

-t^AQ 2 = R(Q,b) (31) 

In addition, AQ* = oAQi + AQ 2 for all a. If more than one trial value of a is needed in the 
one-dimensional search, Eqs. (30) and (3 1 ) clearly provide a more efficient means of finding a 
proper update of AQ* than Eq. (29). 
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3.2.2 Implementation of SAADO 


Computationally, to assist in the implementation of the line-search procedure of 
Kqs. (26}-(29), a merit function is constructed using the interior-penalty-function method. This 
merit function is a combination of the objective, side-constraint, and aerodynamic residual equa- 
tions and is given by 

111 

4>(Q, b, r) = b) + r ^ 1 (Q, b) + r' 1/2 R T (Q, b) R (Q, b) (32) 

where r is the penalty-function parameter (which is initially assumed to be unity and is 
progressively reduced from one design iteration cycle to the next). 1 he value of o, as 
determined by the line-search process, is designed to ensure that the numerical value of the 
merit function is reduced when compared with its value from the previous design-iteration cycle 
(i.e., 4>(CT,bV) < <h(Q,b,r)). 

At the final optimum design, Q* and b* are required to satisfy not only the Kuhn-Tucker 
conditions but must also reduce the residual of the aerodynamic equations to a very small 
value; that is, Q* should represent a high-quality solution of the steady-state aerodynamic field 
equations. One procedure for meeting this condition is to perform some flow-solver iterations to 
reduce the residual of the flow equations somewhat (thus, Q* is improved and b* is held fixed) 
before the next design-optimization iteration is started. This procedure can also be supplemented 
by using an appropriate weighting coefficient tor the aerodynamic field equations in the merit 
function, at each cycle in the design optimization process. These and other issues in regard to 
the implementation of the SAADO algorithm are discussed later in greater detail. 
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If the steady-state flow equations (Eq. (3)) are linear in Q, then a single Newton iteration 
will always solve these equations exactly; in this case, Eqs. (24) and (25) become identical 
to Eqs. (4) and (5), respectively. Furthermore, after Ab is found, Eq. (11) can be solved 
directly for the new Q* without resorting to the approximation process described by Eq. (29). 
Therefore, if the state equations are linear, then the SAADO scheme is exactly equivalent to the 
conventional optimization method and no computational savings is achieved. 

The proposed SAADO procedure can be summarized as follows: 

Step 1. Select an initial design b and assume an initial value of the held variables Q (e.g., 
the free-stream value). 

Step 2. Set up and solve Eqs. (20) and (23) for the adjoint variables A 0 and Aj, i = 1, 
2, . . . , m. 

Step 3. Set up the linearized optimization problem of Eqs. (24) and (25), and solve for Ab 
with a linear programming method. 

Step 4. Solve Eq. (29) or Eqs. (30) and (31) for AQ\ 

Step 5. Find the step size or with the line-search procedure described by Eqs. (26) — (3 1), 
including use of the merit function (Eq. (32)); update the design and field variables. 

Step 6. Check all convergence criteria; if the criteria are not satisfied, then return to Step 
2 and use a smaller value for the penalty parameter r. 
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3.3 Extension to Multidimensional Flow Problems 


The basic SAADO algorithm is quite general and, in principal, can be applied to any flow 
problems. However, the successful extension of the SAADO algorithm to those flow problems 
depends upon the particular CFD solvers used and how they are employed in iteratively solving 
nonlinear flow equations such as Eq. (3), as well as linear equations such as (Eqs. (20), (23), 
and, (29)). For one-dimensional and some two-dimensional flow problems, a Newton direct 
solver can be used to solve these equations on current computers. However, a strict application 
of Newton’s method is not feasible on modem computers when solving these equations for 
large two- and three-dimensional problems with the Euler or Navier-Stokes equations. An exact 
construction of the true Jacobian coefficient matrix and a direct in-core Newton solution 
strategy applied to the nonlinear flow equation (such as Eq. (3)) at each iteration is impossible on 
current supercomputers for these problems because of insufficient memory. This same difficulty 
also applies, of course, to the direct solution of the large linear systems (Eqs. (20), (23) and 
(29)) associated with the proposed SAADO scheme. However, iterative solution of these large 
linear systems of equations is practical for multidimensional problems by recasting them in 
“incremental” or “delta” form. The details of this incremental formulation and a demonstration 


of the feasibility and benefits for aerodynamic problems are given in references [11-13]. 
The incremental iterative form, for example, of Eq. (20) is 



(Ao) m +l — (^o) m + (AA 0 ) m 


(33) 
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where m is the incremental iterative index. An equation similar to Eq. (33) is obtained for 
Eq. (23). The incremental iterative form of Eq. (29) is 

-(i) A(AQ,)l '" - + ®) (nAb) + R(Q ’ b) 

(ACn m+ i = (AQX + [A(AQ-)] m (34) 

Note that in Eqs. (33) and (34) the tilde over the Jacobian operator on the left-hand side 
indicates that this operator can be approximate, whereas those on the right-hand side cannot. 
Most three-dimensional CFD flow solvers typically do not construct the exact Jacobian matrix 
operator the idea here is to use the efficient left-hand-side operator of the existing flow 
solver and construct derivative terms on the right-hand side via automatic differentiation [14, 
15]. Automatic differentiation has been applied to advanced CFD codes [16, 17] to obtain 
consistent, discrete sensitivity derivatives of aerodynamic output functions with respect to input 
and modeling parameters. It has not yet been applied, however, in an incremental iterative form, 
where significant computational gains are expected. 

For the present implementation and demonstration of SAADO, a one-dimensional flow 
problem will be considered. Therefore, a Newton direct-solver method will be used to solve the 
linear adjoint equations, as well as to iterate on the nonlinear flow problems. 

4.0 OPTIMIZATION PROBLEMS AND NUMERICAL RESULTS 

Two quasi-one-dimensional sample problems have been selected to demonstrate the SAADO 
scheme; its performance on these problems will be compared with that of the more conventional 
“black box” optimization scheme. One of the objectives of this study is to ensure that the 
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SAADO algorithm is capable of delivering a high-quality flow solution at the end of the design- 
optimization process. The effect of applying flow-solver type iterations between optimization 
cycles will be studied to achieve this objective. 

The first problem is a fully supersonic flow through a diverging nozzle; the second is a 
transonic flow through a converging-diverging nozzle, with a normal shock in the diverging 
section. These sample problems are described subsequently in more detail. 

4.1 Problem Statement 

In each of the two quasi-one-dimensional sample problems, a “target nozzle” shape is defined 
through the a priori specification of a target value of each of the design variables; then, a fully 
converged, steady-state numerical solution of the governing equations is obtained for the target 
nozzle. Thus, the target flow speed iij is known at each jth grid point of the target-nozzle 
geometry. Then, an initial nozzle shape (different, of course, from that of the target nozzle) 
is defined by specification of an initial value for each of the design variables. In the sample 
problems to be presented, the objective of the optimization procedure is to manipulate the nozzle 
shape by varying the design variables such that an optimum match is obtained between the 
flow-speed profile that evolves during the optimization process and the flow-speed profile of the 
target nozzle. The objective function »P 0 is defined by using the difference between the (current) 
evolving flow speed and the target flow speed at each grid point; that is, 

jmax 

E (Bj - (35) 

1 j=i 

where uj is the flow speed that evolves at the jth grid point and jmax is the total number of grid 
points. No side constraints are present in the sample problems. 
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The flow speed uj is obtained by solving the quasi-one-dimensiona! Euler equations, which 


are 


r)Q tfF(Q) 

ic + fc 1 +S(Q)= 0 


(36) 


where Q = [p, pu, />e 0 ] 1 ,F(Q)= [pu, pu 2 + P, (pe G + P)u] r , and S(Q) = -$£y[pu, 
/>u 2 , (pe 0 -f P)u] T . In these equations, p is density, u is flow speed, P is pressure, e 0 is the 
specific total energy (i.e., e 0 = e + u 2 / 2, where e is the specific internal energy), and A(x) is the 
local cross-sectional area. The ideal-gas law with a constant ratio of specific heats 7 (taken to 
be 1.4) is used for closure, which implies P = (7 — 1 )(pe 0 — pu 2 / 2). The governing equations 
are discretized and solved numerically with the upwind flux-vector-splitting method of Van Leer 
[18], which includes the use of higher order accuracy to approximate the flux terms. A more 
complete discussion of these numerical procedures is presented in reference [19]. The flow 
field was discretized with 100 grid points. For steady flow, the discretization of the governing 
equations, together with the numerical treatment of the boundary conditions, results in a large 
set of coupled nonlinear algebraic equations with the form of Eq. (3). In this study, the basic 
procedure for solving the discrete nonlinear flow equations is Newton’s root-finding method, or 
a modified version of the same, which is described subsequently in this section. 


4.2 Sample Problem 1: Supersonic Nozzle 

The streamwise variation of the cross-sectional area A(x) for the fully supersonic nozzle 
problem is given by 


A(x) = b| + l>2 tanli (0.8 x — !>,-$) (0 < x < 10) (37) 


18 



where b ( , b 2 , and b 3 are selected to be the design variables (i.e., b = (bi , b 2 , b3 ) ). The initial 
values of these design variables are bj = 1.33, b 2 = 0.30, and b 3 = 3.90; their target values are 
b| = 1.398, b 2 = 0.347, and b 3 = 4. For this fully supersonic flow problem, the Mach number on 
the inflow boundary is 1.5; p, u, and P are specified and held fixed (at the free-stream values) on 
this boundary and are extrapolated on the outflow boundary. In constructing the merit function 
(Eq. (32)) for this problem, weighting coefficients were set to 1 (unity) for both the objective 
function and the root-mean-square value of the residual of the aerodynamic equations. The 
penalty-function parameter r (Eq. (32)) was updated at every optimization step by halving it. 

In the computational results that follow, the new optimization scheme is compared with the 
conventional black-box approach, wherein the design optimization module and the flow-analysis 
module are totally independent. The design optimization software used for this was Automated 
Design Synthesis (ADS) [20]. The conventional design optimization approach was employed 
with gradients (derivatives of the objective function with respect to the design variables) 
calculated by both one-sided finite difference (FD) and quasi-analytical (QA) differentiation of the 
flow-analysis code. The method of feasible directions was used as the optimizer, and the golden- 
section method was used for the one-dimensional search for the two cases above that involve the 
conventional design optimization approach. In addition to the basic SAADO algorithm explained 
above, two variants of this algorithm were also applied to the present sample problem to study the 
effect of introducing additional flow-solver iterations between the optimization cycles. The first 
variation deviates from the basic SAADO algorithm in the following manner; several Newton 
iterations (three in this study) are performed on the initial guess for the field variables (which 
in this case is the free stream), and then the actual design-optimization process can be started 
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and continued as explained in section 3.2.2. This process means that the three Newton iterations 
are introduced immediately after step 1 (section 3.2.2). This particular procedure is denoted 
herein as algorithm 1. In the second variation, a few Newton iterations are introduced into 
the basic SAADO algorithm after Q* (the updated field variables) are obtained to improve the 
quality of the solution for b* (the updated design variables) held fixed. This means that a few 
Newton iterations are introduced immediately after step 5 (section 3.2.2). This procedure is 
denoted herein as algorithm 2a if one Newton iteration is added; it is denoted algorithm 2b if 
two Newton iterations are added. 

All results presented here were computed on a Cray-2; the flow-analysis code requires nine 
Newton iterations for a complete, fully converged flow-analysis (steady-state) solution. The 
design optimization algorithms were applied for 20 cycles each. Figure 1 compares the initial- 
and target-velocity profiles for the nozzle (although initial values for field variables at interior 
points are free-stream values). Table 1 compares the algorithms in terms of the final objective 
function, the root-mean-square value of the residual of the aerodynamic equations at the end of 
the design process, the total number of equivalent Newton-Raphson (NR) iterations, the overall 
computational efficiency (i.e., total CPU time), and the final design variables achieved. From 
the results, algorithm 2b clearly gives the best overall results in terms of the magnitude of the 
final objective function and the quality of the final flow solution. The other algorithms also give 
good results because for most engineering purposes, a residual value of 1.0E-4 is considered 
to be satisfactory. These results show that the number of intermediate Newton iterations and 
their placement in the SAADO algorithm enables the user to control the residual and, hence, the 
quality of the flow solution at the end of the design process. 
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To study the performance of the basic algorithm, the residual and the objective function 
after each optimization cycle are plotted in Fig. 2. The detailed result shows that a global 
reduction occurs in the residual as the optimization progresses, but the overall residual at the end 
of the design process is considered satisfactory for engineering purposes. The effectiveness of 
using intermediate Newton iterations (algorithms 2a and 2b) in achieving the previously stated 
objectives is thus reiterated. 

Figure 3 shows the residual history during the optimization process using algorithm 2b 
superimposed on a single standard flow-solver residual history. In this plot, each peak represents 
the beginning of an optimization cycle; each valley represents the end of a cycle. For the basic 
algorithm, Fig. 4 shows the residual change due to changes in the first, second, and third design 
variables; this residual change is substantially different than that in the similar plot shown in 
Fig. 5, which is generated with algorithm 2b. The boxed points in these figures represent the 
residual values at the end of an optimization cycle. A decreasing trend can be seen in the peaks 
in Fig. 5, although it is not clearly distinguishable. 

The final velocity profiles obtained by the SADDO algorithms are all in very good agreement 
with the target-velocity profile. However, note that the sample problem does not impose any 
constraint on the mass flow or the size of the nozzle and, hence, may not yield a unique solution. 
That is, the optimization algorithm (depending on the initial design) can converge to any of the 
local minimum points that are characteristic of the objective function. For example, in Table 1, 
the final values for design variables b are changed in the direction opposite to that which would 
bring them closer to their target values for all methods. Several different initial design variables 
were tested with the same unconstrained objective function and SAADO scheme (algorithm 2b). 
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Each resulted in a different final design; the results are tabulated in Table 2. Figure 6 shows the 
area profiles for all final designs that were generated, as well as the area profile generated for 
the target. In all cases, the final flow-speed profile was extremely close to the target profile. 

4J Sample Problem 2: Transonic Nozzle 

The variation of the cross-sectional area A in the streamwise direction x in this sample 
problem is given by the expression 

A(x) = 1 — b|x + Ihx 2 (0 < x < 1) (38) 

where b| and f >2 are the design variables. The target-design variables are b| = 0.8 and f >2 = 
0.8, and their initial values are 0.9 and 1 .2, respectively. For this transonic flow problem, the 
frce-stream Mach number is 0.55, and a normal shock is found in the diverging part of the 
nozzle. On the inflow boundary, the entropy and stagnation enthalpy are fixed at the free-stream 
values; the flow speed is extrapolated. On the outflow boundary, the density and flow speed 
are extrapolated; the static pressure is held fixed (i.e., P/Poo = 0.9 in the present sample). 
Weighting coefficients of 20 for the objective function and 50 for the residual were set in the 
merit function (Eq. (32)). The penalty-function parameter is updated in a manner similar to 
the first sample problem. 

With a strict application of Newton’s root-finding method to solve the nonlinear flow 
equations, the method commonly “overshoots” the root and diverges. This overestimation is 
particularly true during the initial start-up phase of the solution procedure and is most commonly 
experienced when higher order accurate spatial discretizations are applied to flow fields with 
shocks. Typically, this problem can be overcome with the use of underrelaxation. In this 
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study, an underrelaxation factor was included to initiate the flow-solver iterations; however, this 
initial use of underrelaxation was progressively reduced during the iterative process in proportion 
to the reduction in the average value of the residual of the flow equations. Thus, the use of 
underrelaxation was progressively removed during the solution process, such that as the root (i.e., 
the steady-state flow solution) was approached the strict use of Newton’s root-finding method 
was recovered. 

In this second sample problem, without the use of underrelaxation, the strict implementation 
of the basic SAADO scheme (as outlined in section 3.2) was divergent, which was corrected 
with modifications to the basic algorithm. First, basic Newton flow-solver iterations with 
undcrrelaxation were implemented in a manner as described previously to initially reduce the 
residual of the flow equations approximately two orders of magnitude (1.0E-2). Then, the 
modified SAADO algorithm was begun by evaluating the adjoint-variable vectors) (i.e., by 
solving Eq. (20) (and Eq. (23) if side constraints were present)) directly without the inclusion of 
any relaxation factor. After the solution of the linear programming problem (Eqs. (24) and (25)) 
for Ab, Eq. (29) is solved directly for AQ* without any relaxation factor. Then, two Newton 
flow-solver iterations, including the use of an underrelaxation factor that was progressively 
removed (as described previously), were introduced to reduce the residual immediately after 
the design changes (i.e., immediately after step 4 of section 3.2.2). After the design and state 
variables were updated, three Newton-type flow-solver iterations (including the progressively 
removed underrelaxation factor) were again used to improve the solution, and the new design 
b* was held fixed (i.e., immediately after step 5 of section 3.2.2). 
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Figure 7 compares the initial-, final-, and target-velocity profiles. The SAADO design 
optimization was carried out for 78 cycles, and its results are reported in Table 3, including 
the objective function, residual value, and efficiency of the scheme in terms of CPU time 
and equivalent flow-solver iterations. The residual history and objective function during the 
design process are shown in Fig. 8. From these results, the intermediate flow-solver steps 
appear to provide a useful tool to control the quality of the flow solution at the end of the 
optimization process. Although not explicitly shown here, note that depending on several 
different combinations of weighting coefficients and the number of intermediate flow-solver 
type iterations used the final results can vary slightly in terms of the optimum design and the 
quality of the final flow solution that is achieved. 

In Fig. 9, the residual history during the optimization process is shown superimposed on a 
single standard flow-solver residual history. Figure 10 shows a sample of these residual values 
at each design iteration on an expanded scale. The peaks and the valleys represent the beginning 
and the end, respectively, of each optimization cycle. The decreasing trend in the residual is 
easily discernible in Fig. 10. 

5.0 CONCLUSIONS 

An optimization technique for nonlinear aerodynamic problems, simultaneous aerodynamic 
analysis and design optimization (SAADO), has been investigated in this study. This algorithm 
enables the simultaneous reduction of the objective function and the correction of the constraint 
violations; in this strategy, the constraints include the residuals of the aerodynamic equations. 
This design optimization technique requires that the steady-state aerodynamic equations be 
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satisfied only at the final optimum solution of the design optimization problem. This requirement 
obviates the need to obtain the complete the steady-state flow solution accurately at every 
optimization iteration, which is characteristic of conventional design optimization procedures and 
accounts for the bulk of central processing unite (CPU) time in the entire design optimization 
process. 

The SAADO algorithm and its variants are successfully implemented on two quasi-one- 
dimensional nozzle-flow problems, including a sample problem where a discontinuity exists (i.e., 
a normal shock) in the flow field. A substantial increase occurs in the number of design variables 
and constraints when the field variables are treated as part of the set of design variables and when 
the aerodynamic flow equations are considered to be part of the set of constraints; this is nullified 
by the introduction of the adjoint-variable vector. Furthermore, the appropriate introduction of 
direct Newton iterations (including the use of an underrelaxation factor in the design optimization 
process) helps to control the value of the residual of the aerodynamic equations. Depending upon 
the number of flow-solver iterations used and their placement in the design optimization scheme, 
a trade-off can be achieved between the values of the residual of the aerodynamic equations (at 
the end of the design process) and the objective function and the CPU time. In both sample 
problems, a considerable savings in terms of equivalent Newton iterations (and, hence, overall 
CPU time) is achieved with the SAADO scheme when compared with the conventional black- 
box design optimization procedure. 

In principle, the SAADO algorithm is extendable to large-scale multidimensional flow 
problems. However, in these cases, the incremental iterative method should be employed to 
solve the linear equations that are involved in this new optimization scheme. 
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Table 1. Optimization Results: Supersonic Nozzle Problem 
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Table 2. Illustration of Nonuniqueness of Optimum Design Problem 
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Table 3. Results for Transonic Flow Problem 
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Figure 3. Residual history of SAADO algorithm 2b and 
flow-analysis code (supersonic flow). 
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Figure 4. Trend of residual with changes in design (basic algorithm). 


35 



Or 



[o] Residual values a! 
end of optimization 
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Figure 5. Trend of residual with changes in design (SAADO algorithm 2b). 


36 




Figure 6. Area profiles of different final designs generated by SAADO algorithm 2b 
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Figure 7. Initial-, final- and target- velocity profiles (transonic flow). 
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Figure 8. Residual and objective function (transonic flow). 
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Figure 9. Residual history of SAADO 
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Figure 10. Trend of residual with changes in design (transonic flow). 
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